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Abstract. Case's method is generalized to higher spatial dimensions by reformulat- 
ing the method of rotated reference frames. 



1. Introduction 

The method of solving the radiative transport equation with elementary solutions or 
singular eigenf unctions is called Case's method [1]. Although the method gives insight 
into the theoretical structure of the specific intensity, it works only for one spatial 
dimension [2j. On the other hand, another method called the method of rotated reference 
frames [3l H] has been developed as an efficient numerical algorithm for higher spatial 
dimensions. In this paper, we extend Case's method to three dimensions by considering 
singular eigenf unctions in rotated reference frames. The calculation presented here 
turns out to be another formulation of the method of rotated reference frames. Thus, 
the generalization in this paper is not formal but numerically feasible. 

The reminder of this paper is organized as follows. In section [21 we introduce the 
radiative transport equation. In section [3l we obtain normal modes. In section HI we 
consider eigenvalues. Section |5] is devoted to the method of rotated reference frames. 
In section El we obtain the three-dimensional Green's function in an infinite medium. 
In section [71 we present the Green's function for isotropic scattering. Finally, we give 
summary in section [SI Polar and azimuthal angles in rotated reference frames are 



computed in Appendix A The expansion coefficients in the method of rotated reference 



frames are calculated in [Appendix B 



2. Radiative transport equation 

Let u{r,s) be the specific intensity at position r G M^ in direction s G S^. The (time- 
independent) radiative transport equation is given by 

s ■ VM(r, s) + ^(r, s) = c / f{s-s')u{r,s')ds'. (1) 

We assume < c < 1. Let /i = cos6' be the cosine of the polar angle of s and (p be the 
azimuthal angle of s. Suppose the scattering phase function /(s ■ §') can be modeled by 
a polynomial of spherical harmonics of order A^. 

N I 
/=0 m=-l 

Let g = /g2(s ■ §')/(§ ■ §') ds' and /g2 /(§ ■ §') ds' = 1. We have N = 0, fo = 1 in the 
case of isotropic scattering, N = 1, fo = 1, fi = g in the case of linear scattering, and 
N = CO, fi = g'' in the case of the Henyey-Greenstein model [5]. Let us express /(s • s') 
as 

1=0 m=-l ^ '' 

Here the polynomials p™(/i) are related to PJ^ifi) as [6] 



They satisfy the following recurrence relations and orthogonality relations. 

(/ - m + l)prti(/^) = (2/ + l)/ipr(/^) -(1 + m)pT.ii^^): (5) 

1 2(l-\-mV 

where we introduced 

dm(/i) = (l-^')'""' d/i. (7) 

3. Normal modes 

We seek solutions of the form of plane-wave decomposition [7J [S]. We introduce i/ G M 

and q G M^, and define vector k G C'^ as 

k = -k, k = I '""^ ] , Q(Hq|) = v^r + ^^W- (8) 

We emphasize that k and k are functions of u and q. We assume the specific intensity 
of the form 

«-(r, s; q) = $-(§; k)e-^-^ = ^i^, /i(k)) (l - ^(k)^) '"'^' e'-^«e-'^•^ (9) 

where //(k) and fCk) are the cosine of the polar angle of s and the azimuthal angle of 
s, respectively, in the rotated reference frame whose ;2-axis coincides with the direction 



of k (See Appendix A). We normalize 0*" as 



— y^0"^(i/,/i(k))(l-Mk)') ds = y 0™(i.,/i)dm(/i) = l. (10) 

By plugging the ansatz (^ into the radiative transport equation ([T]), we obtain 



Mk)K^,.. ,„^^^^ ,.nN2\l"^l/2 



imip(k) 



z/,/i(k))(l-Mkf) e 

c^^ / (s(k) . s'(k)) 0'"(z/,/i'(k)) (l - /i'(k)2)'"''^%^-^'(^) ds', (11) 

where we used s ■ k = /i(k) and expressed /(§ ■ s') in the rotated reference frame. The 
right-hand side is calculated as 



11/2 

RHS = 2ttcQ {N - \m\) (l - i^{k)A "" e''^^^'^) ^ /; 



l-l/^ir..(fc) ^ , 2r + l{l'-m) 
'' 47r (/' + m)! 



1 



X pI7(/i(k))y ^p-(/.')0™(^,/^')dm(/i'), (12) 

where the step function B(-) is defined as Q{x) = 1 for x > and = for x < 0. Hence, 

N 



(z.-Mk))r(^,Mk)) = 27rcz.e(Ar-|m|) J^ f, J |^^ ^ ^j, p[r(/i(k))/^[^(^), (13) 



where we defined 

hT{v) = y^^0-(z/,/x)pr(/i)dm(/i). (14) 

Since the right-hand side of (iT3l) is zero for |r?7,| > A^ and then (j)"^ = 0, hereafter we 
suppose 

< |m| < N. (15) 

Let us define 

ai = l-cfiQ{N-l). (16) 

From (iT3l) . we obtain 

aiuhriiy) = J /i0"^(i/,/i)pr(/^) dm(/i). (17) 

Thus the recurrence relation for h^^i^) is obtained [9]: 

u{2l + l)aihf{u) - (/ - m + l)Kll^{u) - {I + m)Kl'_^{u) = 0, (18) 

with h'^^iu) = (2m)!/2™m! for m > and h'^^iu) = (-l)™/2l™l(|m|!) for m < 0, and 
hl'^^^^iu) = (2|m| + l)(z/ + a|„| - l)/i|"|(z/). We also have 

/,-H(^)^(_1)h|[^M|/,H(^). (19) 

The functions hY'iy) are computed using (fT8l) . 
Let us define 

7'"(z/,/i(k)) = J2 (2^' + l)/r§^p[^(/i(k))/.ir(z.). (20) 



l'=\m\ 



(/' + m) 



We note that 7 *"(//, /i(k)) = 7'"(z/, /i(k)). The function c/)'" is obtained as 

r{y. /i(k)) = ^p2:!>lii^ + A-(z.)(l - z.^)-H<5(z. - ^(k)). (21) 

^ z/ — /i(k) 

4. Discrete eigenvalues and continuous spectrum 

By integrating fl^Tl) with respect to s, we obtain 

l = ^v! 2!^i^dm(/i)+ /" A"^(z/)5(z/-/i)d/i. (22) 

2 j_i u — n J_i 

Let us define 

A-(^) = l-f r^^^^^^^dm(/i), (23) 

where 2; G C. Eigenvalues z/ ^ [~1) 1] are solutions to 

A™(z/) = 0. (24) 



We write these discrete eigenvalues as i^^j" (t'o" > J-'i' > ■ • ■ > ^Ti-i > !)• Note that 
z/""^ = u"^. The number of discrete eigenvalues M depends on |m| and we have fiU\ |6] 

M < N - \m\ + I. 

For u G [—1, 1], we have 

X"i.) = l-'-^Vr'-^i,nM. (25) 

Note that A-™(z/) = A™(z/). This implies 0-™(z/,/i(k)) = 0™(z/,/i(k)). 

In two dimensions (r G M^, § G S"*^), we have s = (cosy^, siny)) and /(§ ■ s') = 
^ + :^ X]m=i /"i'^°^^('^ ~ V^')- We can introduce elementary solutions as ^(r, s) = 
0(z/, /i(k))e"^''", where /i(k) = cos(v9 — y^g.). Here k = k/z/ and k = (cosy^j^, siny^j.) = 
{—\vq,Q{yqj) with g G M. Let us define hmiy) = /o'^cos(m(y9)0(z/, /i) dy? and 

7(1/, /i(k)) = X]m=-Ar /m cos fmy9(k) j /iml'^)- The singular eigenfunctions in two 
dimensions are obtained as 

0(z., ^^{k)) = |lp2i^^^ + A(z.)5(z. - /.(k)), (26) 

where X{v) = 1 — f^ ^/o'' ^^^-^ d*/?. Discrete eigenvalues are obtained as solutions to 

A(^) = i_^ r"2l£ii^d<^ = o, zee. (27) 

27r Jo 2; - /i 
In particular, for isotropic scattering A^ = 0, we have h.{z) = 1 — |f /o'^~~dv9 = 
1 — T^T = 0- Therefore, we obtain z = ±l/-\/l — c^ = iz/Q. 

5. Method of rotated reference frames 

Let us expand singular eigenfunctions with spherical harmonics. By introducing c™(z/), 

we write 

00 

$-(s;k)= 5^cr(z/)l^„.(s;k). (28) 

l=\m\ 

The calculation of the specific intensity by this expansion is called the method of rotated 
reference frames [3lll]. From flTTj) . we obtain 

cr(^) --Y.il f^Yi'mi^Li^) ds) q?(//) = cfiQiN - l)cT{v). (29) 

Hence we arrive at an eigenproblem: 

B^\^"^{u)) = u\^^{u)), (30) 

where 

verier I, y§2 



p _ ^72^ / (/ + 1)2 — m^ 

-Si',i-i + W ..., , ^.„ TT (^r,i+i, (31) 



(4/2 - l)a,(T,_i -'-^ V (4(^ + 1)' - l)^/+i^; 



m-{u)) 



-^icT{v) 
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(32) 



where the normahzation Z™(z/) will be determined below so that (■?/'™'(z/)| ■?/'"'(//)) = 1 is 
satisfied. Note that v depends on \m\. We obtain 



$ir(s;k) = yz^ E ^^^^y..(s;k). 



l=\m\ 



(33) 



6. The Green's function 



Let us consider the Green's function of the radiative transport equation in an infinite 
medium. The Green's function obeys 

s ■ VG'(r, s; ro, Sq) + G(r, s; fq, Sq) = c / /(s ■ s')G(r, s'; fq, Sq) ds' 

+ 5(r-ro)<5(s-So). (34) 

To proceed, we introduce $™(s; k) such that 



§2 



/i$™(s;k) $:;)(§; k 



ds = byijl, 



(35) 



where the Kronecker delta b^yi is understood as the Dirac delta b{y — u') for the 
continuous spectrum. The function $™(s;k) will be determined as we compute the 
Green's function. 

We replace the source term in (134|) by a jump condition and solve 



§ • VG'(r, s; ro, Sq) + ^(r, s; fq, Sq) = c / /(s ■ s')G'(r, s'; fq, Sq) ds', 

J §2 



(36) 



G{p, Zq + 0, s; ro, Sq) - G{p, Zq - 0, s; fq, Sq) 



s ■ z 



'^(P - Po)^ (s - So) 



with the boundary condition lim|r|_>oo G'(r, s; tq, Sq) = and r = (p, z) (p G M^, 5; G 
Let us expand the Green's function using normal modes ([9]). 

N „ rM-1 



G{p,z,s;Pq,Zo,So) = V / y'af+(q)Mj"V(r,s;q 

™_ ^r 7m2 L ._n 



A-(qX(r,s;q)dz/ 



m=-Af ""^ '-j=0 

dq 



N „ rA/-l 



Z > Zq, 



G{p, z, s; po, zo, §o) = - 5^ / E af_(q)M™ (r, s; q 



(37) 



A-(q)C(r,s;q)dz/ 



m=-N ""^ L j=0 

dq 



(2vr) 



2' 



Z < Zo- 



From the jump condition, coefficients a"i and A^ are determined as 






A'"(q) = Q-'^'i'PoQQi^l'iD^o/i' 



*r±(§o;k) 
$r(§o;k) 



(38) 
(39) 



Hence the Green's function is written as 

±1 



CIp, 2;,s;po,zo,so) 



(2vr) 



^iq-(p-Po) 



TV 

E 

m=-N 



M-1 



5^<l>fi(s;k) l>7±(so;k) 



g-QK^IqDk-^ol/^T" 



j=o 



where upper signs are chosen for z > zq and lower signs are chosen for z < zq- 

Noting the relation fl5^ . the Green's function obtained with the method of rotated 
reference frames |1] is expressed as 



N 



j<i-{p-Po) 



EE 



1 



X $±,(s;k) [$![?,(so;k)l*e-^("l''l)l"~"°l/"dq. 



G'(r,s;ro,So) 



(41) 
By comparing ( HOl) and ( 14T]) . we obtain 

<l^,(s; k) = $™,(s; k) [±z.g(Hq|)^'"(/^)]-^ ■ (42) 

To determine Z™(z/), we consider the one- dimensional case. By integrate the 
Green's function over Pq, we obtain 



G{Z, S; Zq, S 



.) = ±E{ 





m=-Af '^ 


M-1 


E^Ji(§'^) 


<^r±(§o;z) 


i=o 





-k-^olK 



+ $?.(§; z 



^^.(so;^)] e-l^-^«l/^dz/ . (43) 



On the other hand, the one- dimensional Green's function is given by [TOl |6] 



1 ^ r 

G (2;, s; 2:0, So) = ^ X] j 

'Kn _ \T ^ 



m=-N 



M-1 



E TLr(±<,/i)0'"(±^r,/io) (1 -/i2)'"^'e-i^-^«i/''" 



..0 ^. 



^^' xa^'^"^^''' '")'^'"(±^' -"o) (^ - -"') '"' ^"''"'"'^'^ 



du 



where 



1 c\\"^(7) 



J 2''^ ' "'^ ''^ ' dz 

^^{p) = z/A™+(z/)A™-(z/) (1 - z/^)"'""' . 
Here A™'^(i/) = lime_s>o+ A™(z/ ± ie). By comparing (I43p and flU]) . we obtain 
i/™Z'"(i/f ) = 27rAr™, z/Z"(z/) = 27rA/'(z/). 



(44) 

(45) 
(46) 

(47) 



Finally we obtain the Green's function as 

1 



1 



(2vr) = 



^iq-(p-Po) 



TV 

E 

m=-N 



M-1 



i=o "^^"J 



g(z/f I q|)A/;" 



1 



*r±(§o;k) 



^-Q{i,--\ci\)\z-zo\/u^ 



$-(s;k) 



Qi^hlWO 
7. Isotropic scattering 

For the isotropic scattering A^ = 0, we have 

<^°(§; k) = ^P -^ + A°(z.)5(z. - /x(k)) = <p{u, /x(k)), 

where 



(48) 



(49) 



d/i = 1 - CI/ tanh-i z/. (50) 

.1 u — fi 



Here /i(k) is given in Appendix A In this case M = 1 and the discrete eigenvalues 

,0 




±i/q = ±z/o are solutions to 



AO(^) = 1 - ^ 



1 1 

d/x = 1 — cztanh"^ - = 0. 

2 y_i z — fj, z 



(51) 



The Green's function is obtained as 



G{p,z,s;P(),Zo,So) 



{2ny 



MiP-Po) 



^.±z/o,Mk)U* ±^o,/io(k) e-«('^°l''l)l^-^''l/'^° 
g(z/o|q|)A/o V / V J 

(±u, /i(k)) 0* (±u, /io(k)) e-«("l'il)l"-"°l/'^ dz/| dq. 



+ 



1 



Ql^lqDA/"!^) 






(52) 



8. Summary 



We have constructed a three-dimensional version of Case's method of solving the 
radiative transport equation. Extending Case's method to higher spatial dimensions was 
considered in 1960's shortly after Case's paper [1]. However the generalization is only 
formal and cannot be used for numerical calculation [21[TT1[7]. In this paper, the singular 
eigenfunctions are generalized using rotated reference frames and they are labeled by 
Case's discrete eigenvalues and continuous spectrum. The method developed in this 
paper is equivalent to the method of rotated reference frames. Therefore, the present 
method can also be considered a numerical method. Indeed, the three-dimensional 
radiative transport equation with anisotropic scattering in an infinite medium [3], H] and 



9 

in the slab geometry [12] were numerically calculated. The radiative transport equation 
is also solved in two dimensions |13l [m [151 iMl [H]- The use of the method of rotated 
reference frames or the Green's function (HHI) for inverse problems is also proposed |18j . 

Appendix A. Polar and azimuthal angles in rotated reference frames 

Let 6 and (p be the polar and azimuthal angles of s in the laboratory frame. Let yjg. and ^j^ 
be the polar and azimuthal angles of k in the laboratory frame. For k = (— iz/q, (5(z^|q|)), 
we obtain 

cos^j^ = k ■ z = (5(z^|q|), sin6'5. = yl^^cos^^ = i|z/q|, (A.l) 

and 

{(Pa + vr z/ > 0, 

(A.2) 
(y9q Z/ < 0, 

where y^q is the angle of q. Therefore, we have 

;u(k) = s ■ k = — iz^|q| sin 6 cos((y9 — (y9q) + Ql'^lql) cos 0. (A. 3) 

In general, we can rotate functions as follows. Let us introduce rotated spherical 
harmonics Yjm(s;k) [3]: 

Yi^is- k) = V{y)Yi„,{s) = Y, Dl,^iip^, 9^, 0)F,„,(s), (A.4) 

m'=—l 

where Dl^,^{ip^^,6^^,0) = e~''^''^^d'-^,^{6^). Here rf^,^ are the Wigner d-matrices [19]. 
That is, Yim{s,h) are spherical harmonics defined in a rotated reference frame whose 
z-axis coincides with the direction of the unit vector k. We have Yim{s) = Yim{s] z). We 
write analytically continued Wigner's (i-matrices as 

c^!n'm(^k) = dL'm[^'r{i'h\)]- (A.5) 

First a few matrices are obtained as 

We note that d^^, = (-l)-+-'ci'_^_^, = (- 1) '"+'"' c?^,^- All c?^,^[ir(z/|q|)] are 
computed using the recurrence relations [I2]. We obtain 



e-.M -(^- „„:,,V"" H-irVM2m+l) 



(l-Mfc)-) '"^Cm S »-"--CJ''.)>w(s). (A.7) 



(2m + 1)!! 
where ^ satisfies cos 6 = fi with /x in 



m'=—m 
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Appendix B. Expansion coefficients 

Here we calculate cj'iv). We have 

2\|m|/2 imip 



2\-|"i| 



+ A™(z/)(l-z/2)-'™'5(z.-/i) 



(l-/.2)i'"i/^e-^F,:„(s)ds. 



(B.l) 



Hence, 



crfz/) = 27r, 



'2/ + 1 (/-m)! 
47r (/ + m)! 



TE/'"(2r + i)j;;^W 



r' = |m| 



(Z" + m)! 



/I 
5(z/ - /i) d/i 



(B.2) 



Note that cj"(-z/) = (-l)'+™c["(z/) because Pr{-y) = (-l)'+'"P;'"(z/). Therefore, we 
obtain for u ^ [—1,1] 

cT{v) = 2vri 



/2/ + l(/-^)!c. ^^^^^^2^„^^^(r-^^„^^^ 



47r (/ + m)! 2 ^ 



{I" + m)\' '" 
X (-l)'"2g™,,(,,,,)(z/)P-„(,,,,)(z/), (B.3) 

where <5max(n")('^) ^^"-^ -^min(n")('^) have a branch cut from — oo to 1 [20], and for 

i^G[-l,l] 

TV 



crfz/) = 27r 



'2/ + 1 {l-m)\ 



'-^i:f'^'i-"^^^'^;i^-. 



l"=\m\ 



{I" + m)\ 



z/ — e 



id 



+ A'"(z/)(l-z/2)-H/2p^'^ 



(B.4) 
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